# Load Phyloseq object generated in HIV_oral_preprocessing.Rmd
# Options include GreenGenes, SILVA, RDP and HitDB
ps0 <- readRDS("~/Dropbox/Research/human_volunteer/data/ps0.Harris.rdp.RDS")
ps0
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1936 taxa and 122 samples ]
## tax_table() Taxonomy Table: [ 1936 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 1936 tips and 1934 internal nodes ]
# Load mapping file
map <- import_qiime_sample_data("~/Dropbox/Research/human_volunteer/data/mapping_human_volunteer.txt")
ps0 <- merge_phyloseq(ps0, map)
ps0
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1936 taxa and 122 samples ]
## sample_data() Sample Data: [ 122 samples by 119 sample variables ]
## tax_table() Taxonomy Table: [ 1936 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 1936 tips and 1934 internal nodes ]
# View sample variables & generate basic stats
sample_variables(ps0)
## [1] "X.SampleID" "BarcodeSequence"
## [3] "LinkerPrimerSequence" "sequence_ID"
## [5] "miseq_run" "plate"
## [7] "well" "collaborator"
## [9] "tissue" "number_seqs_slout"
## [11] "sample_surviving_dada2" "sample_ID"
## [13] "tube_ID" "patient_ID"
## [15] "randomization" "day"
## [17] "randomization_arm" "Date_assessment"
## [19] "Age" "Gender"
## [21] "Ethnicity" "Asthma"
## [23] "Eczema" "Hayfever"
## [25] "Allergies" "Def_allergies"
## [27] "Other_illness" "Chronic_medication"
## [29] "Def_medication" "Herbs_OTC"
## [31] "def_Herbs_OTC" "Alcohol"
## [33] "Alcoholunitsweek" "Smoking"
## [35] "Drugs" "Drugstype"
## [37] "Drugsusage" "Weight"
## [39] "Length" "BMIcalc"
## [41] "SystolicBP" "DiastolicBP"
## [43] "Pulse" "Respfreq"
## [45] "Physical_exam_abnorm" "def_phys_abnormalities"
## [47] "Lab_abnorm_screen" "def_abnorm_lab_screen"
## [49] "Serum_screen" "PBMC_screen"
## [51] "Feces_screen" "Date_day_9"
## [53] "New_medication_day_9" "def_new_med_day_9"
## [55] "Course_completed" "Antibiotic_side_effects"
## [57] "Def_abx_side_effects" "Lab_abnorm_post_antibiotics"
## [59] "def_abnorm_post_antibiotics" "Date_day_0"
## [61] "Feces_day_0" "New_medication_day0"
## [63] "Def_new_medication_day_0" "Vaccination_side_effects"
## [65] "Def_vac_side_effects" "Lab_abnorm_post_vaccination"
## [67] "def_abnorm_post_vaccination" "Date_day_7"
## [69] "Feces_day_1" "Feces_day_2"
## [71] "Feces_day_3" "Feces_day_4"
## [73] "Feces_day_5" "Feces_day_6"
## [75] "Feces_day_7" "Serum_day_7"
## [77] "New_medication_day7" "def_new_med_day7"
## [79] "Date_day_14" "New_medication_day14"
## [81] "def_new_med_day14" "Serum_day_14"
## [83] "PBMC_day_14" "Date_day_28"
## [85] "New_medication_day28" "def_new_med_day28"
## [87] "Serum_day_28" "PBMC_day_28"
## [89] "Study_complete" "pre_RV_IgA"
## [91] "d7_RV_IgA" "boost_change"
## [93] "d7_RV_IgA_boost_fc" "d7_roto_boost"
## [95] "d7_rota_boost_updated" "d14_RV_IgA"
## [97] "d28_RV_IgA" "RV_Ag_shedding_d0"
## [99] "RV_Ag_shedding_d1" "RV_Ag_shedding_d2"
## [101] "RV_Ag_shedding_d3" "RV_Ag_shedding_d4"
## [103] "RV_Ag_shedding_d5" "RV_Ag_shedding_d6"
## [105] "Shedding" "pneumo_IgG_pre"
## [107] "pneumo_IgG_7" "pneumo_IgG_14"
## [109] "pneumo_IgG_28" "pneumo_IgG_7vPre"
## [111] "pneumo_IgG_14vPre" "pneumo_IgG_28vPre"
## [113] "tetanus_IgG_pre" "tetanus_IgG_d7"
## [115] "tetanus_IgG_d14" "tetanus_IgG_d28"
## [117] "tetanus_IgG_7vPre" "tetanus_IgG_14vPre"
## [119] "tetanus_IgG_28vPre"
sd(sample_sums(ps0))
## [1] 23974.15
get_taxa_unique(ps0, "Phylum")
## [1] "Bacteroidetes" "Firmicutes"
## [3] "Verrucomicrobia" "Proteobacteria"
## [5] "Fusobacteria" "Synergistetes"
## [7] "Actinobacteria" "Spirochaetes"
## [9] NA "Lentisphaerae"
## [11] "Euryarchaeota" "Elusimicrobia"
## [13] "Cyanobacteria/Chloroplast" "Tenericutes"
# Create a new data frame of the sorted row sums, a column of sorted values from 1 to the total number of individuals/counts for each RSV and a categorical variable stating these are all RSVs.
readsumsdf = data.frame(nreads = sort(taxa_sums(ps0), TRUE),
sorted = 1:ntaxa(ps0),
type = "ASVs")
# Add a column of sample sums (total number of individuals per sample)
readsumsdf = rbind(readsumsdf,
data.frame(nreads = sort(sample_sums(ps0), TRUE),
sorted = 1:nsamples(ps0),
type = "Samples"))
# Make a data frame with a column for the read counts of each sample for histogram production
sample_sum_df <- data.frame(sum = sample_sums(ps0))
# Make plots
# Generates a bar plot with # of reads (y-axis) for each taxa. Sorted from most to least abundant
# Generates a second bar plot with # of reads (y-axis) per sample. Sorted from most to least
p.reads = ggplot(readsumsdf, aes(x = sorted, y = nreads)) +
geom_bar(stat = "identity") +
ggtitle("ASV Assessment") +
scale_y_log10() +
facet_wrap(~type, scales = "free") +
ylab("# of Reads")
# Histogram of the number of Samples (y-axis) at various read depths
p.reads.hist <- ggplot(sample_sum_df, aes(x = sum)) +
geom_histogram(color = "black", fill = "firebrick3", binwidth = 150) +
ggtitle("Distribution of sample sequencing depth") +
xlab("Counts") +
ylab("# of Samples")
# Final plot, side-by-side
grid.arrange(p.reads, p.reads.hist, ncol = 2)
# Basic summary statistics
summary(sample_sums(ps0))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 1205 20670 24463 41152 99585
Interpretation: Lots (20+) samples with zero ASV’s. These will likely need to be removed in the subsequent step. Other than that the distributions look normal. Skew from high to low ASV is due to the incorporation of antibiotics treated samples which result in fewer ASV.
# Format a data table to combine sample summary data with sample variable data
ss <- sample_sums(ps0)
sd <- as.data.frame(sample_data(ps0)) # useful to coerce the phyloseq object into an R data frame
ss.df <- merge(sd, data.frame("ASV" = ss), by ="row.names") # merge ss with sd by row names. Rename ss to ASVs in the new data frame
# Plot the data by the treatment variable
y = 100 # Set a threshold for the minimum number of acceptable reads. Can start as a guess
x = "randomization_arm" # Set the x-axis variable you want to examine
label = "Description" # This is the label you want to overlay on the points that are below threshold y. Should be something sample specific
# Plot
p.ss.boxplot<- ggplot(ss.df, aes_string(x, y = "ASV")) + # x is what you assigned it above
geom_boxplot(outlier.colour="NA", aes(group = randomization_arm)) +
scale_y_log10() +
geom_hline(yintercept = y, lty = 2) + # Draws a dashed line across the threshold you set above as y
geom_jitter(alpha = 0.6, width = 0.15, size = 2.5) +
#geom_text(data = subset(ss.df, RSVs <= y), aes_string(x, y="RSVs", label=label), size=2, hjust -1) +# This labels a subset that fall below threshold variable y and labels them with the label variable
ggtitle("Number of ASV") +
ylab("ASV (log10)") +
facet_grid(~day) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(legend.position = "NULL") +
theme(axis.title.x = element_blank())
p.ss.boxplot
# Remove samples with fewer than 100 ASV
# Save samples that fell below the selected threshold
ps.failed <- prune_samples(sample_sums(ps0) < y, ps0)
nsamples(ps.failed)
## [1] 22
# Create a table of failed sample information
failed.samples <- data.frame(sample_names(ps.failed), sample_data(ps.failed)$sample_ID, sample_data(ps.failed)$tube_ID)
write.table(failed.samples, file = "./results/failed_samples.txt", sep = "\t")
# Remove samples with fewer than y ASV
ps0
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1936 taxa and 122 samples ]
## sample_data() Sample Data: [ 122 samples by 119 sample variables ]
## tax_table() Taxonomy Table: [ 1936 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 1936 tips and 1934 internal nodes ]
ps0 <- prune_samples(sample_sums(ps0) >=y, ps0)
ps0
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1936 taxa and 100 samples ]
## sample_data() Sample Data: [ 100 samples by 119 sample variables ]
## tax_table() Taxonomy Table: [ 1936 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 1936 tips and 1934 internal nodes ]
min(sample_sums(ps0))
## [1] 103
Interpretation: There does not appear to be a relationship between low numbers of ASVs, treatment or time-point. So removing samples that behave unexpectedly (< 100 sequences).
# Outlier evaluation
out.bray <- ordinate(ps0, method = "NMDS", distance = "bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1817636
## Run 1 stress 0.1867117
## Run 2 stress 0.1883909
## Run 3 stress 0.1836108
## Run 4 stress 0.1833127
## Run 5 stress 0.1889399
## Run 6 stress 0.1870656
## Run 7 stress 0.1925691
## Run 8 stress 0.1880074
## Run 9 stress 0.1927114
## Run 10 stress 0.1895174
## Run 11 stress 0.188155
## Run 12 stress 0.1821197
## ... Procrustes: rmse 0.02361735 max resid 0.1442131
## Run 13 stress 0.1965863
## Run 14 stress 0.18175
## ... New best solution
## ... Procrustes: rmse 0.000991212 max resid 0.009441048
## ... Similar to previous best
## Run 15 stress 0.1919634
## Run 16 stress 0.1966386
## Run 17 stress 0.184672
## Run 18 stress 0.1870022
## Run 19 stress 0.2003686
## Run 20 stress 0.1821406
## ... Procrustes: rmse 0.02418486 max resid 0.1442341
## *** Solution reached
p.NMDS.outlier <- plot_ordination(ps0, out.bray, color = "randomization_arm", axes = c(1,2)) +
theme_bw() +
facet_grid(~day) +
geom_point(size = 2) +
ggtitle("NMDS of Bray Distances \nOutlier Evaluation") +
stat_ellipse(type = "norm", geom = "polygon", alpha = 1/10, aes(fill = randomization_arm))
p.NMDS.outlier
Interpretation: There is no distiguishable feature to identify samples as obvious outliers using this projection. No samples will be removed based on NMDS of Bray-curtis distances.
# Begin by removing sequences that were not classified as Bacteria or were classified as either mitochondria or chlorplast
ps0 # Check the number of taxa prior to removal
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1936 taxa and 100 samples ]
## sample_data() Sample Data: [ 100 samples by 119 sample variables ]
## tax_table() Taxonomy Table: [ 1936 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 1936 tips and 1934 internal nodes ]
ps1 <- ps0 %>%
subset_taxa(
Family != "mitochondria" &
Class != "Chloroplast"
)
ps1 # Confirm that the taxa were removed
## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 1590 taxa and 100 samples ]
## sample_data() Sample Data: [ 100 samples by 119 sample variables ]
## tax_table() Taxonomy Table: [ 1590 taxa by 7 taxonomic ranks ]
## phy_tree() Phylogenetic Tree: [ 1590 tips and 1588 internal nodes ]
saveRDS(ps1, file = "./results/ps1.RDS")
# Intermediate data subsetting for prevelance plotting
# Full subsetting with removal of zero count data will be performed in a subsequent chunk
ps1.d_9 <- subset_samples(ps1, day == "-9")
ps1.d0 <- subset_samples(ps1, day == "0")
ps1.d7 <- subset_samples(ps1, day == "7")
# Prevelance estimation
# Calculate feature prevelance across the data set
prevdf.d_9 <- apply(X = otu_table(ps1.d_9),MARGIN = ifelse(taxa_are_rows(ps1.d_9), yes = 1, no = 2),FUN = function(x){sum(x > 0)})
prevdf.d0 <- apply(X = otu_table(ps1.d0),MARGIN = ifelse(taxa_are_rows(ps1.d0), yes = 1, no = 2),FUN = function(x){sum(x > 0)})
prevdf.d7 <- apply(X = otu_table(ps1.d7),MARGIN = ifelse(taxa_are_rows(ps1.d7), yes = 1, no = 2),FUN = function(x){sum(x > 0)})
# Add taxonomy and total read counts to prevdf
prevdf.d_9 <- data.frame(Prevalence = prevdf.d_9, TotalAbundance = taxa_sums(ps1.d_9), tax_table(ps1.d_9))
prevdf.d0 <- data.frame(Prevalence = prevdf.d0, TotalAbundance = taxa_sums(ps1.d0), tax_table(ps1.d0))
prevdf.d7 <- data.frame(Prevalence = prevdf.d7, TotalAbundance = taxa_sums(ps1.d7), tax_table(ps1.d7))
# Create a table of Phylum, their mean abundances across all samples, and the number of samples they were detected in
plyr::ddply(prevdf.d_9, "Phylum", function(df1){cbind(mean(df1$Prevalence),sum(df1$Prevalence))})
## Phylum 1 2
## 1 Actinobacteria 2.000000 92
## 2 Bacteroidetes 2.308594 591
## 3 Elusimicrobia 0.500000 1
## 4 Euryarchaeota 3.000000 6
## 5 Firmicutes 2.067068 2404
## 6 Fusobacteria 0.200000 3
## 7 Lentisphaerae 1.400000 7
## 8 Proteobacteria 1.282609 118
## 9 Spirochaetes 0.000000 0
## 10 Synergistetes 0.000000 0
## 11 Tenericutes 0.000000 0
## 12 Verrucomicrobia 2.800000 14
plyr::ddply(prevdf.d0, "Phylum", function(df1){cbind(mean(df1$Prevalence),sum(df1$Prevalence))})
## Phylum 1 2
## 1 Actinobacteria 1.8478261 85
## 2 Bacteroidetes 1.2265625 314
## 3 Elusimicrobia 0.0000000 0
## 4 Euryarchaeota 3.0000000 6
## 5 Firmicutes 1.1160791 1298
## 6 Fusobacteria 0.4666667 7
## 7 Lentisphaerae 0.8000000 4
## 8 Proteobacteria 2.0543478 189
## 9 Spirochaetes 1.0000000 1
## 10 Synergistetes 1.5000000 3
## 11 Tenericutes 1.0000000 1
## 12 Verrucomicrobia 3.2000000 16
plyr::ddply(prevdf.d7, "Phylum", function(df1){cbind(mean(df1$Prevalence),sum(df1$Prevalence))})
## Phylum 1 2
## 1 Actinobacteria 2.347826 108
## 2 Bacteroidetes 2.007812 514
## 3 Elusimicrobia 0.500000 1
## 4 Euryarchaeota 3.000000 6
## 5 Firmicutes 2.480653 2885
## 6 Fusobacteria 1.133333 17
## 7 Lentisphaerae 1.400000 7
## 8 Proteobacteria 2.239130 206
## 9 Spirochaetes 1.000000 1
## 10 Synergistetes 1.500000 3
## 11 Tenericutes 0.000000 0
## 12 Verrucomicrobia 4.400000 22
# Reorder based on visual inspection of relative dominance of each taxa for plotting
prevdf.d_9$Phylum <- factor(prevdf.d_9$Phylum, levels = c("Firmicutes", "Bacteroidetes", "Proteobacteria", "Actinobacteria", "Verrucomicrobia", "Lentisphaerae", "Fusobacteria", "Elusimicrobia", "Spirochaetes", "Synergistetes", "Tenericutes"))
prevdf.d0$Phylum <- factor(prevdf.d0$Phylum, levels = c("Firmicutes", "Bacteroidetes", "Proteobacteria", "Actinobacteria", "Verrucomicrobia", "Lentisphaerae", "Fusobacteria", "Elusimicrobia", "Spirochaetes", "Synergistetes", "Tenericutes"))
prevdf.d7$Phylum <- factor(prevdf.d7$Phylum, levels = c("Firmicutes", "Bacteroidetes", "Proteobacteria", "Actinobacteria", "Verrucomicrobia", "Lentisphaerae", "Fusobacteria", "Elusimicrobia", "Spirochaetes", "Synergistetes", "Tenericutes"))
#Prevalence plot - t1
prevdf1.d_9 <- subset(prevdf.d_9, Phylum %in% get_taxa_unique(ps1.d_9, "Phylum"))
p.prevdf1.d_9 <- ggplot(prevdf1.d_9, aes(TotalAbundance, Prevalence / nsamples(ps1.d_9),color=Family)) +
# geom_hline(yintercept = 0.03, alpha = 0.5, linetype = 2) +
geom_point(size = 1.5, alpha = 0.7) +
scale_x_log10(limits = c(1,10000), breaks = c(100,10000)) +
xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
facet_grid(~Phylum) +
theme(legend.position="None") +
theme(axis.title.x = element_blank()) +
ggtitle("Day -9")
#Prevalence plot - t2
prevdf1.d0 <- subset(prevdf.d0, Phylum %in% get_taxa_unique(ps1.d0, "Phylum"))
p.prevdf1.d0 <- ggplot(prevdf1.d0, aes(TotalAbundance, Prevalence / nsamples(ps1.d0),color=Family)) +
# geom_hline(yintercept = 0.03, alpha = 0.5, linetype = 2) +
geom_point(size = 1.5, alpha = 0.7) +
scale_x_log10(limits = c(1,10000), breaks = c(100,10000)) +
xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
facet_grid(~Phylum) +
theme(legend.position="None") +
theme(axis.title.x = element_blank()) +
ggtitle("Day 0")
#Prevalence plot - t3
prevdf1.d7 <- subset(prevdf.d7, Phylum %in% get_taxa_unique(ps1.d7, "Phylum"))
p.prevdf1.d7 <- ggplot(prevdf1.d7, aes(TotalAbundance, Prevalence / nsamples(ps1.d7),color=Family)) +
# geom_hline(yintercept = 0.03, alpha = 0.5, linetype = 2) +
geom_point(size = 1.5, alpha = 0.7) +
scale_x_log10(limits = c(1,10000), breaks = c(100,10000)) +
xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
facet_grid(~Phylum) +
theme(legend.position="None") +
theme(axis.title.x = element_blank()) +
ggtitle("Day 7")
grid.arrange(p.prevdf1.d_9, p.prevdf1.d0, p.prevdf1.d7, nrow = 3)
Interpretation: There are a number of low prevalent bacteria. There is a notable reduction in both Firmicutes and Bacteroidetes following antibiotic treatment. This is not subset by treatment (randomization_arm), but the overall reduction suggests the antibiotics reduced the prevalence of many taxa following treatment. There is a visible restoration by day 7.
Interpretation: Log10 transformation does tend to shift the data to normality, but it is not striking and may ‘overshift’ the data. Should consider using log10 transformed data for future analysis, but consider other transformations as needed.
# Community composition - Baseline and Week 24 timepoints combined
# Phyla level plots
# Melt to long format (for ggploting)
# Prune out phyla below % in each sample
# Creat a Phyla level table for external analysis
# write.table(ps1_phylum, file = "./results/phylum.txt", sep="\t")
threshold = 0.03
ps1_phylum <- ps1 %>%
tax_glom(taxrank = "Phylum") %>% # agglomerate at phylum level
transform_sample_counts(function(x) {x/sum(x)} ) %>% # Transform to rel. abundance
psmelt() %>% # Melt to long format
filter(Abundance > threshold) %>% # Filter out low abundance taxa
arrange(Phylum) # Sort data frame alphabetically by phylum
# Convert Sample No to a factor because R is weird sometime
ps1_phylum$patient_ID <- as.factor(ps1_phylum$patient_ID)
p.comm.bar <- ggplot(ps1_phylum, aes(x = patient_ID, y = Abundance, fill = Phylum)) +
geom_bar(stat = "identity", width = 0.9) +
facet_wrap(randomization_arm~day, scales = "free_x") +
ggtitle("Community Composition Over Time and Randomization Arm") +
ylab("Relative Abundance") +
theme(legend.position = "bottom") +
theme(axis.text.x = element_blank()) +
theme(axis.title.x = element_blank()) +
scale_fill_brewer(palette = "Dark2")
p.comm.bar
# Draw interactive plot
ggplotly(p.comm.bar)
Interpretation: Similar to the prevalence plots, there is a notable alteration in the Bacteroidetes and Firmicutes in antibiotic treated samples at Day 0. This is still noticeable at Day 7, but the Firmicutes are making somewhat of a comeback. Proteobacteria invade following antibiotic treatment as well, but are largely absent, particularly in the narrow spectrum antibiotic samples after 7 days. Overall, the restoration to the pre-antibiotic bacterial microbiome appears to be more prominent following narrow spectrum treatment than in broad spectrum treatment.
my.comparisons.paired <- list(c("Pre", "0"), c("Pre", "7"), c("0","7"))
p.paired.rich <- ggpaired(ps1.rich, x = "day", y = "richness_0", outlier.shape = NA, id = "patient_ID") +
geom_jitter(width = 0.2) +
ylab("Richness") +
theme(axis.title.x = element_blank()) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
facet_grid(~randomization_arm) +
stat_compare_means(label = "p.signif", method = "t.test", ref.group = "-9", hide.ns = TRUE) +
theme(axis.text.x = element_blank())
p.paired.sd <- ggpaired(ps1.rich, x = "day", y = "diversities_shannon", outlier.shape = NA, id = "patient_ID") +
geom_jitter(width = 0.2) +
ylab("Shannon Diversity") +
theme(axis.title.x = element_blank()) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
facet_grid(~randomization_arm) +
stat_compare_means(label = "p.signif", method = "t.test", ref.group = "-9", hide.ns = TRUE)
grid.arrange(p.paired.rich, p.paired.sd, nrow = 2)
Interpretation: Richness (the number of taxa) and Shannon diversity (relative abundance and evenness) are similar in all untreated sampels over time. Both narrow and broad spectrum antibiotic treatment decrease alpha diversity over time, with some recovery occurring in both treatement protocols between Day 0 and 7.
IgA boost defined as a greater than or equal to 2 fold change between pre-rotavirus IgA levels and Day 7 IgA levels.
p.rich.rota_boost <- ggboxplot(ps1.rich, x = "d7_rota_boost_updated", y = "richness_0", outlier.shape = NA) +
geom_jitter(width = 0.2) +
ylim(0,250) +
ylab("Richness") +
theme(axis.title.x = element_blank()) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
stat_compare_means(label = "p.signif", method = "t.test") +
facet_grid(randomization_arm~day) +
ggtitle("Rota-IgA Boost") +
theme(axis.text.x = element_blank())
p.sd.rota_boost <- ggboxplot(ps1.rich, x = "d7_rota_boost_updated", y = "diversities_shannon", outlier.shape = NA) +
geom_jitter(width = 0.2) +
ylim(0,5) +
ylab("Shannon diversity") +
theme(axis.title.x = element_blank()) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
stat_compare_means(label = "p.signif", method = "t.test") +
facet_grid(randomization_arm~day)
grid.arrange(p.rich.rota_boost, p.sd.rota_boost, nrow = 2)
Interpretation: There are signficiantly more types of bacteria (richness) at day 7 in both narrow spectrum and broad spectrum antibiotic patients. There is increased bacterial diversity in day -9 pre-treatment samples in the patients that go on to not receive antibiotics, but the levels are normalized over the course of the study.
Shedding was defined as having one or more stool samples per patient positive for rotavirus shedding.
p.rich.shedding <- ggboxplot(ps1.rich, x = "Shedding", y = "richness_0", outlier.shape = NA) +
geom_jitter(width = 0.2) +
ylim(0,250) +
ylab("Richness") +
theme(axis.title.x = element_blank()) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
stat_compare_means(label = "p.signif", method = "t.test") +
facet_grid(randomization_arm~day) +
ggtitle("Shedding") +
theme(axis.text.x = element_blank())
p.sd.shedding <- ggboxplot(ps1.rich, x = "Shedding", y = "diversities_shannon", outlier.shape = NA) +
geom_jitter(width = 0.2) +
ylim(0,5) +
ylab("Shannon diversity") +
theme(axis.title.x = element_blank()) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
stat_compare_means(label = "p.signif", method = "t.test") +
facet_grid(randomization_arm~day)
grid.arrange(p.rich.shedding, p.sd.shedding, nrow = 2)
Interpretation: Both richenss and diversity were increased in day 0 samples of patients treated with narrow spectrum antibiotics.
# Test for age dependency on alpha diversity
# Age per randomization_arm
ancova(richness_0~Age*randomization_arm, data = ps1.rich)
## Analysis of Variance Table
##
## Response: richness_0
## Df Sum Sq Mean Sq F value Pr(>F)
## Age 1 87 87.3 0.0309 0.860816
## randomization_arm 2 28660 14329.9 5.0724 0.008092 **
## Age:randomization_arm 2 3549 1774.6 0.6282 0.535799
## Residuals 94 265557 2825.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ancova(diversities_shannon~Age*randomization_arm, data = ps1.rich)
## Analysis of Variance Table
##
## Response: diversities_shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## Age 1 0.809 0.8093 1.1815 0.2798223
## randomization_arm 2 13.361 6.6806 9.7530 0.0001416 ***
## Age:randomization_arm 2 0.296 0.1480 0.2160 0.8061022
## Residuals 94 64.388 0.6850
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# BMI
ancova(richness_0~BMIcalc*randomization_arm, data = ps1.rich)
## Analysis of Variance Table
##
## Response: richness_0
## Df Sum Sq Mean Sq F value Pr(>F)
## BMIcalc 1 2372 2372.3 0.8461 0.36000
## randomization_arm 2 27040 13520.0 4.8222 0.01015 *
## BMIcalc:randomization_arm 2 4891 2445.3 0.8722 0.42140
## Residuals 94 263550 2803.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ancova(diversities_shannon~BMIcalc*randomization_arm, data = ps1.rich)
## Analysis of Variance Table
##
## Response: diversities_shannon
## Df Sum Sq Mean Sq F value Pr(>F)
## BMIcalc 1 0.372 0.3717 0.5545 0.4583477
## randomization_arm 2 13.049 6.5243 9.7334 0.0001439 ***
## BMIcalc:randomization_arm 2 2.425 1.2126 1.8090 0.1694713
## Residuals 94 63.009 0.6703
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Alpha diversity is dependent on randomization arm, but not BMI
Interpretation: Alpha diversity is not dependent on age or BMI.
## All sample - Treatment
# Richness
p.gam.rich.treatment <- ggplot(ps1.rich, aes(x = day, y = richness_0, color=d7_rota_boost_updated)) +
geom_point(alpha = 0.6, size = 3) +
stat_smooth(method = "gam", formula = y ~ s(x, bs = "cr", k = 3)) +
ggtitle("General Additive Model (GAM) - Treatment") +
facet_grid(~randomization_arm) +
ylab("Richness")
# Shannon Diversity
p.gam.sd.treatment <- ggplot(ps1.rich, aes(x = day, y = richness_0, color=d7_rota_boost_updated)) +
geom_point(alpha = 0.6, size = 3) +
stat_smooth(method = "gam", formula = y ~ s(x, bs = "cr", k = 3)) +
ggtitle("General Additive Model (GAM) - Treatment") +
facet_grid(~randomization_arm) +
ylab("Shannon Diversity")
grid.arrange(p.gam.rich.treatment, p.gam.sd.treatment, nrow = 2)
## All sample - Treatment
# Richness
p.gam.rich.shedding <- ggplot(ps1.rich, aes(x = day, y = richness_0, color = Shedding)) +
geom_point(alpha = 0.6, size = 3) +
stat_smooth(method = "gam", formula = y ~ s(x, bs = "cr", k = 3)) +
ggtitle("General Additive Model (GAM) - Treatment") +
facet_grid(~randomization_arm) +
ylab("Richness")
# Shannon Diversity
p.gam.sd.shedding <- ggplot(ps1.rich, aes(x = day, y = richness_0, color=Shedding)) +
geom_point(alpha = 0.6, size = 3) +
stat_smooth(method = "gam", formula = y ~ s(x, bs = "cr", k = 3)) +
ggtitle("General Additive Model (GAM) - Treatment") +
facet_grid(~randomization_arm) +
ylab("Shannon Diversity")
grid.arrange(p.gam.rich.shedding, p.gam.sd.shedding, nrow = 2)
# Ordination
# weighted unifrac
ord.ps1.wuni <- ordinate(ps1, method = "PCoA", distance = "wunifrac")
ord.ps1.wuni.log <- ordinate(ps1.log, method = "PCoA", distance = "wunifrac")
# unifrac
ord.ps1.uni <- ordinate(ps1, method = "PCoA", distance = "unifrac")
# Bray-curtis
ord.ps1.bray <- ordinate(ps1, method = "PCoA", distance = "bray")
# PCoA plot
p.ord.wuni <- plot_ordination(ps1, ord.ps1.wuni, color = "randomization_arm") +
geom_point(size=3, alpha = 0.7) +
scale_color_brewer(palette = "Dark2") +
geom_point(colour = "grey50", size = 1.5) +
facet_grid(~ day) +
stat_ellipse(geom = "polygon", alpha = 1/15, aes(fill = randomization_arm)) +
labs(color = "randomization_arm") +
ggtitle("wUniFrac")
p.ord.wuni.log <- plot_ordination(ps1.log, ord.ps1.wuni.log, color = "randomization_arm") +
geom_point(size=3, alpha = 0.7) +
scale_color_brewer(palette = "Dark2") +
geom_point(colour = "grey50", size = 1.5) +
facet_grid(~ day) +
stat_ellipse(geom = "polygon", alpha = 1/15, aes(fill = randomization_arm)) +
labs(color = "randomization_arm") +
ggtitle("wUniFrac: log transformed")
p.ord.uni <- plot_ordination(ps1, ord.ps1.uni, color = "randomization_arm") +
geom_point(size=3, alpha = 0.7) +
scale_color_brewer(palette = "Dark2") +
geom_point(colour = "grey50", size = 1.5) +
facet_grid(~ day) +
stat_ellipse(geom = "polygon", alpha = 1/15, aes(fill = randomization_arm)) +
labs(color = "randomization_arm") +
ggtitle("UniFrac")
p.ord.bray <- plot_ordination(ps1, ord.ps1.bray, color = "randomization_arm") +
geom_point(size=3, alpha = 0.7) +
scale_color_brewer(palette = "Dark2") +
geom_point(colour = "grey50", size = 1.5) +
facet_grid(~ day) +
stat_ellipse(geom = "polygon", alpha = 1/15, aes(fill = randomization_arm)) +
labs(color = "randomization_arm") +
ggtitle("Bray-Curtis")
grid.arrange(p.ord.wuni, p.ord.wuni.log, p.ord.uni, p.ord.bray, nrow = 4)
Interpretation: Selecting wUniFrac as it 1) tends to agree with alpha diversity observations and 2) explains a high-percentage of variance on axis 1. However, there are likely details to be resolved from each measure.
Community composition differences beetween groups is primarily explained by treatment (randomization arm). As with alpha diversity, the extent of difference between groups is greatest at day 0 and diminished by day 7. However, the differences between each sample at day 7 are larger than at day 7 (more variant across Axis 1 and 2) suggesting intra-individual spread in bacterial communities.
p.ord.boost <- plot_ordination(ps1, ord.ps1.bray, color = "d7_rota_boost_updated") +
geom_point(size=3, alpha = 0.7) +
scale_color_brewer(palette = "Dark2") +
geom_point(colour = "grey50", size = 1.5) +
facet_grid(randomization_arm ~ day) +
stat_ellipse(geom = "polygon", alpha = 1/15, aes(fill = d7_rota_boost_updated)) +
ggtitle("wUniFrac Rotavirus IgA Boost")
p.ord.shedding <- plot_ordination(ps1, ord.ps1.bray, color = "Shedding") +
geom_point(size=3, alpha = 0.7) +
scale_color_brewer(palette = "Dark2") +
geom_point(colour = "grey50", size = 1.5) +
facet_grid(randomization_arm ~ day) +
stat_ellipse(geom = "polygon", alpha = 1/15, aes(fill = Shedding)) +
ggtitle("wUniFrac Rotavirus Shedding")
grid.arrange(p.ord.boost, p.ord.shedding, nrow = 2)
Interpretation: Visual inspection makes it appear that bacterial community structure is different between samples from patients that boosted and those that did not boost at Day 7. There is poor sample distribution of boosters and shedders to analyze no antibiotic controls and too few boosted samples to analyze broad spectrum samples.
# Set a random seed so that exact results can be reproduced
set.seed(1000)
# Function to run adonis test on a physeq object and a variable from metadata
doadonis <- function(physeq, category) {
bdist <- phyloseq::distance(physeq, "wunifrac")
col <- as(sample_data(physeq), "data.frame")[ ,category]
# Adonis test
adonis.bdist <- adonis(bdist ~ col)
print("Adonis results:")
print(adonis.bdist)
# Homogeneity of dispersion test
betatax = betadisper(bdist,col)
p = permutest(betatax)
print("Betadisper results:")
print(p$tab)
}
# Two factor adonis
doadonis2 <- function(physeq, category, category2) {
bdist <- phyloseq::distance(physeq, "wunifrac")
col <- as(sample_data(physeq), "data.frame")[ ,category]
col2 <- as(sample_data(physeq), "data.frame")[ ,category2]
# Adonis test
adonis.bdist <- adonis(bdist ~ col * col2)
print("Adonis results:")
print(adonis.bdist)
}
# Runs Permanov and Homogeneity of dispersion test
# See ?betadisper for more info
# Randomization arm at each time point
doadonis2(ps1, "day", "randomization_arm") # day = **, randomization_arm = ***, interaction = n.s.
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col * col2)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 1 0.5808 0.58085 3.9030 0.03398 0.005 **
## col2 2 2.0554 1.02770 6.9055 0.12025 0.001 ***
## col:col2 2 0.4665 0.23323 1.5672 0.02729 0.106
## Residuals 94 13.9894 0.14882 0.81847
## Total 99 17.0921 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
doadonis(ps1.t1, "randomization_arm") # n.s.
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 2 0.06800 0.033998 1.1857 0.10146 0.316
## Residuals 21 0.60215 0.028674 0.89854
## Total 23 0.67014 1.00000
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.02441588 0.012207942 2.903593 999 0.083
## Residuals 21 0.08829296 0.004204427 NA NA NA
doadonis(ps1.t2, "randomization_arm") # ***
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 2 0.78594 0.39297 16.867 0.52112 0.001 ***
## Residuals 31 0.72223 0.02330 0.47888
## Total 33 1.50817 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.001956101 0.0009780507 0.1497809 999 0.853
## Residuals 31 0.202426115 0.0065298747 NA NA NA
doadonis(ps1.t3, "randomization_arm") # **
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 2 0.6298 0.31489 2.5571 0.11593 0.019 *
## Residuals 39 4.8025 0.12314 0.88407
## Total 41 5.4323 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.06811062 0.03405531 2.5677 999 0.088
## Residuals 39 0.51725553 0.01326296 NA NA NA
# Randomization arm over time
doadonis(ps1.con, "day") # n.s.
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 1 0.06192 0.061919 1.0496 0.03176 0.366
## Residuals 32 1.88783 0.058995 0.96824
## Total 33 1.94974 1.00000
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.02073363 0.010366814 1.086088 999 0.345
## Residuals 31 0.29589795 0.009545095 NA NA NA
doadonis(ps1.narrow, "day") # n.s. (p < 0.1) # likley not-significant due to Day -9 and Day 7 similarity
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 1 0.0931 0.093058 0.90936 0.02605 0.427
## Residuals 34 3.4793 0.102333 0.97395
## Total 35 3.5724 1.00000
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.1248340 0.062417018 9.693034 999 0.001
## Residuals 33 0.2124992 0.006439368 NA NA NA
doadonis(ps1.broad, "day") # ***
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 1 0.2855 0.28545 2.0451 0.06807 0.087 .
## Residuals 28 3.9082 0.13958 0.93193
## Total 29 4.1936 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 2 0.4770339 0.23851694 20.64645 999 0.001
## Residuals 27 0.3119159 0.01155244 NA NA NA
# Interactions between day and rotavirus boost in separate treatment groups
doadonis2(ps1.narrow, "day", "d7_rota_boost_updated") # n.s.
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col * col2)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 1 0.1429 0.14286 0.83695 0.02379 0.500
## col2 1 0.2141 0.21414 1.25454 0.03566 0.253
## col:col2 1 0.1854 0.18545 1.08646 0.03088 0.318
## Residuals 32 5.4621 0.17069 0.90966
## Total 35 6.0045 1.00000
# Rotavirus sheeding in each radnomization arm over time
doadonis2(ps1.narrow, "day", "Shedding") # n.s.
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col * col2)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 1 0.1110 0.110998 0.68401 0.02035 0.603
## col2 1 0.0678 0.067804 0.41783 0.01243 0.814
## col:col2 1 0.0820 0.081964 0.50509 0.01503 0.752
## Residuals 32 5.1928 0.162276 0.95218
## Total 35 5.4536 1.00000
doadonis2(ps1.broad, "day", "Shedding") # *** for day, n.s. for shedding (samples are different over time, but not based on rotavirus shedding)
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col * col2)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 1 0.35058 0.35058 4.2481 0.13581 0.003 **
## col2 1 0.03771 0.03771 0.4569 0.01461 0.834
## col:col2 1 0.04753 0.04753 0.5760 0.01841 0.760
## Residuals 26 2.14568 0.08253 0.83117
## Total 29 2.58150 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Interaction between treatment and rotavirus vaccine boost at each time point
doadonis2(ps1.t1, "randomization_arm", "d7_rota_boost_updated") # n.s.
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col * col2)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 2 0.17375 0.086873 1.20147 0.10114 0.295
## col2 1 0.06167 0.061674 0.85295 0.03590 0.422
## col:col2 1 0.10859 0.108588 1.50178 0.06321 0.184
## Residuals 19 1.37382 0.072306 0.79974
## Total 23 1.71783 1.00000
doadonis2(ps1.t2, "randomization_arm", "d7_rota_boost_updated") # n.s. (samples are differnet by randomization arm only)
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col * col2)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 2 3.5305 1.76524 16.9073 0.52575 0.001 ***
## col2 1 0.0562 0.05624 0.5387 0.00838 0.720
## col:col2 1 0.1006 0.10061 0.9637 0.01498 0.447
## Residuals 29 3.0278 0.10441 0.45089
## Total 33 6.7152 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
doadonis2(ps1.t3, "randomization_arm", "d7_rota_boost_updated") # sasmples are differnt by randomization arm (**) and by rotavirus boost status, but there is no interaction
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col * col2)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 2 0.26196 0.13098 2.4310 0.09178 0.037 *
## col2 1 0.58868 0.58868 10.9262 0.20626 0.002 **
## col:col2 2 0.06386 0.03193 0.5926 0.02238 0.589
## Residuals 36 1.93960 0.05388 0.67958
## Total 41 2.85409 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Double check day 7 for narrow as it looks separated in the plot
ps1.t3.narrow <- subset_samples(ps1.t3, randomization_arm == "Narrow spectrum antibiotics")
sample_data(ps1.t3.narrow)$randomization_arm
## [1] Narrow spectrum antibiotics Narrow spectrum antibiotics
## [3] Narrow spectrum antibiotics Narrow spectrum antibiotics
## [5] Narrow spectrum antibiotics Narrow spectrum antibiotics
## [7] Narrow spectrum antibiotics Narrow spectrum antibiotics
## [9] Narrow spectrum antibiotics Narrow spectrum antibiotics
## [11] Narrow spectrum antibiotics Narrow spectrum antibiotics
## [13] Narrow spectrum antibiotics Narrow spectrum antibiotics
## [15] Narrow spectrum antibiotics
## Levels: Narrow spectrum antibiotics
sample_data(ps1.t3.narrow)$day
## [1] 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7
doadonis(ps1.t3.narrow, "d7_rota_boost_updated") # 0.08. < 0.1 is near significance
## [1] "Adonis results:"
##
## Call:
## adonis(formula = bdist ~ col)
##
## Permutation: free
## Number of permutations: 999
##
## Terms added sequentially (first to last)
##
## Df SumsOfSqs MeanSqs F.Model R2 Pr(>F)
## col 1 0.3965 0.39650 2.6446 0.16904 0.087 .
## Residuals 13 1.9491 0.14993 0.83096
## Total 14 2.3456 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## [1] "Betadisper results:"
## Df Sum Sq Mean Sq F N.Perm Pr(>F)
## Groups 1 0.01071119 0.01071119 0.307049 999 0.571
## Residuals 13 0.45349577 0.03488429 NA NA NA
Interpretation:
Treatment groups are the same at day -9, but signficantly different at days 0 (p = 0.001) and 7 (p = 0.02)
Control and narrow spectrum antibiotic groups do not change signficantly over time (p = 0.356 and 0.486). Broad spectum antibiotics do change over time (p = 0.002). The reason narrow spectrum antibiotics is not significant here is because day -9 and day 7 are almost identical. The low variant day 0 data are likely being washed out.
Independent of randomization arm, bacterial communities were the same regardless of rotavirus boost at days -9 (p = 0.637), 0 (p = 0.597). They were different at day 7 between samples from patients that boosted and those that did not (p = 0.004)
There is no signficiant difference in bacterial community structure when comparing patients who boosted to those that did not boost over time in any randomziation arm (No antibiotics p-value: 0.681; Narrow spectum antibiotics pivalue: 0.404, Broad spectrum antibiotics: p-value: 0.615)
There is no signficiant difference in bacterial community structure when comparing patients who shed rotavirus to those that did not shed rotavirus over time in any randomziation arm (No antibiotics p-value: 0.996; Narrow spectum antibiotics pivalue: 0.684, Broad spectrum antibiotics: p-value: 0.780)